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pLn Abstract 

An alternating direction implicit (ADI) orthogonal spline collocation (OSC) method is de- 
scribed for the approximate solution of a class of nonlinear reaction-diffusion systems. Its efficacy 
is demonstrated on the solution of well-known examples of such systems, specifically the Brusse- 
lator, Gray-Scott, Gicrcr-Mcinhardt and Schnakcnbcrg models, and comparisons arc made with 
other numerical techniques considered in the literature. The new ADI method is based on an 
extrapolated Crank-Nicolson OSC method and is algebraically linear. It is efficient, requiring 
^ at each time level only 0{Af) operations where Af is the mimbcr of unknowns. Moreover, it is 

^ shown to produce approximations which are of optimal global accuracy in various norms, and 

I— —I to possess superconvergence properties. 
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(N| 1 Introduction 

^ — I 

^ In this paper, we consider two-component reaction-diffusion problems of the form: 

^ (1.1) ^-PAu = f(u), {x,y,t)enx{0,T], 

(1.2) 9^=^' i^,y,t) edUx iO,T], 

(1.3) uix,y,0)=g\x,y), ix,y) gU, 

where A = d'^/dx'^ + d'^jdy^ (the Laplace operator), J7 is the rectangle (a\^h\) x (02,62) with 
boundary $7 = 17 U 9^2, u = \u\. %L2^ . V = diag[Z?i, _D2]"^ is the diagonal matrix of diffusion 
constants, f = [/i,/2]"^ where the functions /i and /2 specify the reaction kinetics of the system, 
and d/dn is the outward normal derivative on the boundary dVL. Problems of this type model 
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phenomena arising in various fields such as chemistry, physics, ecology and biology. In particular, 
Turing [37] proposed such a reaction-diffusion problem to explain biological pattern formation. See 
plj for a recent discussion of this pioneering work. Over the last several decades, reaction-diffusion 
models have evolved in different fields and have come to be known by names such as the Brusselator 
model |29], the Gray-Scott model [El [18], the Gierer-Meinhardt model [E], and the Schnakenberg 
model |35j . The derivation of these models is discussed in [8], for example. 

Alternating direction implicit (ADI) methods were first introduced in the context of finite 
difference methods by Peaceman and Rachford [30] over fifty years ago, and continue to be studied 
extensively for the solution of a variety of problems; see, for example, El [TUl [121 EOl 1211 E21 
HOj Wl\ 133] and references in these papers. There has been much work done on the development, 
analysis and implementation of ADI orthogonal spline collocation (OSC) methods for the solution 
of various scalar transient problems, an overview of which is given in jT3]. The purpose of this 
paper is to present an ADI extrapolated Crank-Nicolson (ECN) OSC method (called simply the 
ADI method in the sequel) for the numerical solution of ( |1.1| )-(1.3). This ADI method is based on 



an ADI ECN OSC scheme formulated in [1] for general nonlinear scalar parabolic problems. Some 
of the salient features of the ADI method are the following. 

• It is a second-order accurate in time perturbation of a Crank-Nicolson OSC scheme and 
preserves the accuracy of that scheme. 



By using extrapolation to linearize ( 1.1 ), it is algebraically linear; that is, the algebraic systems 



to be solved at each step of the method are linear. 

Like all ADI methods, it reduces the solution of the multidimensional problem to the solution 
of sets of one-dimensional problems in each coordinate direction at each time step. With stan- 
dard choices of bases for the spline space in the OSC spatial discretization, these problems 
involve almost block diagonal (ABD) linear systems [2l[T3], which can be solved efficiently 



using existing software [11| I12j. Furthermore, since the differential operator in (1.1) has con- 
stant coefficients, the ABD coefficient matrices in each coordinate direction are independent 
of time and can be decomposed only once for the entire problem. 

• Unlike finite element Galerkin methods, OSC does not involve the approximation of integrals. 
Moreover, without post-processing, it provides superconvergent approximations at the nodes 
of the partition of the spatial domain Q on which the spline space is defined. Unlike finite 
difference methods, it yields continuous approximations to the solution and its first derivatives 
of high-order accuracy throughout 0,. 

A brief outline of the paper is as follows. Definitions and basic notation are given in section 
2, followed by a description of the ADI scheme in section 3. In section 4, the results of numerical 
experiments with test problems from the literature involving the Brusselator, Gray- Scott, Gierer- 
Meinhardt and Schnakenberg models are presented, and comparisons are made with other published 
solution techniques. Concluding remarks are given in section 5. 

2 Preliminaries 

Let {xi}^jQ and {yjl^o partitions (in general, non-uniform) of [ai,&i] and [02,^2], respectively, 
such that 

ai = xo < xi < • • • < X <x =bi, 02 = yo < yi < • • • < yjv„-i < Vn^ = ^2- 
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Let A4x and A4y he the spaces of piecewise polynomials of degree < r, r > 3, defined by 

Mx = {vGC\ai,bi]:v\[x^_^,x^]GPr, l<i<Nx}, 

My = {vG C7i [02,62] : ^^1 G ^r, 1 < J < Ny}, 

where Pr denotes the set of polynomials of degree < r. Let A4 = Mx ^My, the set of all functions 
that are finite linear combinations of products of the form where £ A4x and S 

Aiy Let {Afc}^~^ be the nodes of the (r — l)-point Gauss quadrature rule on [0, 1] and define 



n - /ta; T, N^,r-1 q _ ft?/ t, "y.''-^ where 
— i^i,fc/i=i,fc=i> yy — t?,;/,=i.i=i) wnere 



Ny,r-l 



ilk = xi_i + /if Afc, = + ^)^h hf 



Xi 



h) = yj-yj-i- 



Then 



is the set of collocation points in il. Let {t„}^j^g be a partition of [0, T] such that tn = nr, where 
r = T/Nt, and, for n = 1, . . . , A'^t — l, let in+i/2 = (^n + in+i)/2. Throughout this paper, for any 
function of t, = <j){-,tn). 



3 The ADI Extrapolated Crank-Nicolson OSC Method 

In this section, we describe the ADI scheme. As we shall see, at each time level, the approximate 
solution is a piecewise polynomial of degree < r in one variable only along lines through the 
collocation points. At the final step (or at any user-specified intermediate step), the approximation 
is converted to a two dimensional piecewise polynomial on Q. 

We compute the OSC approximation = [uh,i,Uh,2]'^ G ^Ax^A to the solution u of ( 1.1 )-( 1.3 ) 
at the final time T in the following way. 



Step la. First, we calculate starting values by interpolation at Gauss points. 

• For each G Qx, we choose u^(.^^, •) G A4y x A4y (along vertical lines through the collocation 
points) such that 

(3.1) uO(r,f^) = g°(r,f^), ^yegyu{a2,b2}. 

• For each G Gy, we choose u^(-, C^) G Mx X-Mx (along horizontal lines through the collocation 
points) such that 

(3.2) Q^(e,a = g°(e,e^), eGg.u{ai,6i}. 



Step lb. We also require a second-order accurate approximation to the exact solution u(x, y, ti/2)) 
which we define in the following way. Using Taylor's theorem, (1.1) with t = and ( |1.3| ), we have, 

u(x,y,ti/2) = g°(x,y) + ^[f(g0(x,2/))+PAg0(x,y)] + O(r2), {x,y) G Q. 

'—1/2 

Then, for each G Gx, we choose (^^, •) G A4y x A4y (along vertical lines through the colloca- 
tion points) such that 

T 



(3.3) 



5i/'(e,e^) = g°(e,a + 



f(uO(r,e^))+p 



5y2 



e (^G 



J" 



5u 



1/2 



dy 



^(r,a) = 0, a = 02,62. 
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With standard choices of bases for M.^ and My, the interpolation processes (3.1), (3.2) and 
(3.3) involve ABD linear systems of orders Ny{r — l)+2, Nx{r — l)+2 and Ny{r — l) + 2, respectively, 
each of which can be solved at a cost of 0{Nx) or 0{Ny); see 

Next, we advance in time; for each n = 0,1, Nt — 1, we perform Step 2 followed by Step 
3 as follows: 

Step 2. For each G Gy, determine u"+^/^(-, ^2^) G x such that 



u 



n+l/2 
h 



r/2 



V 



' n2 n+l/2 



da 



+ 



o n+l/2 

dx 



where S;;+^/^(e, ■)€MyxMy is given by 



3u^(r,-)-ur^(r,-)]/2, n = l,...,Ar,-l, 



(3.4) 

and u^(^^,-), uj^''^(^^,-) are given by (3.1), (3.3), respectively. 

This step comprises a set of independent OSC two-point boundary value problems (TPBVPs) along 
horizontal lines through the collocation points Gy, each of which gives rise to an ABD linear system 
of order Ny{r - 1) + 2. 



Step 3. For each G Gx, determine u^"'"-^(^^, •) G My x My such that 



u 



„+l n+l/2 



h 



U 



h 



r/2 



V 



Qr^2 Qy2 



h (^^■^a) = o^ a = 02, 62, 
(9y 



with as in (|3^). 

This step comprises a set of independent OSC TPBVPs along vertical lines through the collocation 
points Gx, each of which gives rise to an ABD linear system of order Nx{r — 1) + 2. 



Once these two steps are completed, we proceed to the final phase of the algorithm, the pur- 
pose of which is to convert the one dimensional approximations along vertical lines through the 
collocation points Gx, •) G My x My, G Gx, to the desired two dimensional approximate 

solution u/i(-,-,r) G X 7W on J7. To this end, we observe that u/i(^^,-,T) G My x My for 
G Gx, and, for fixed G U {02,62}, we introduce u/i(-,^^,r) G Mx x Mx- Then, to obtain 
the desired approximate solution, we perform the following steps. 



Step 4a. For each G U {02,^2}, determine Uh{',i^ ,T) G Mx x Mx along horizontal lines 
through Gy U {a2, 62} such that 



. dx 



{a,S,y,T)=0, a = ai,bi. 



4 



Step 4b. Next, determine u/j(a,-,T) E A4y x Aiy, where a = ai,6i, along vertical lines through 
{ai, 61} such that 

Uh{a, t\T) = \ih{a, e^ T), ^ eGyU {aa, 62}. 

Step 4c. Lastly, determine U/j(-,-,T) E x A^, along horizontal lines through y G [02)^2] such 
that 

(3.5) u,(r,y,T) = uf (r,y), re a., 



where the boundary values for (3.5) are obtained in Step 4b 



Steps 4a and 4c involve the solution of ABD systems of order Nx{r — 1) + 2 whereas Step 4b 
requires the solution of two ABD systems of order Ny{r — 1) + 2. While we have assumed that the 
approximate solution is required only at the final value T, Steps 4a, b, c can be performed at any 
intermediate time level, tn, ^ < n < Nt — 1. 

It can be shown that the total cost of computing the approximate solution u/j using the ADI 
scheme is 0{r'^NxNyNt), where r is the degree of the piecewise polynomials and Nx, Ny and Nf are 
the numbers of subdivisions in the x, y and t directions, respectively. Details of the implementation 
of the method are similar to those of [H Section 2.2] and are omitted. 

4 Numerical Experiments 



In this section, we consider commonly occurring reaction-diffusion models of the form (1.1)-(1.3) 
and compare results obtained by applying the ADI scheme with results appearing in the literature. 
We also demonstrate the global accuracy and superconvergence properties of the ADI scheme. In all 
of the test problems, we consider uniform partitions in the x and y directions with Nx = Ny = N. 

4.1 Brusselator Model 

The system of partial differential equations associated with the Brusselator model [29] comprises 



( l.lD with 

f (u) = [B + u\u2 - (vl + \)ux,Aux - u\u2^ , 

where ^4, B are constants. Various numerical methods have been proposed for its solution. Twizell 
et al., |38j formulated an algebraically linear finite difference method which is second-order accu- 
rate in space and time. Adomian [1] proposed a decomposition method which was subsequently 
corrected and extended in [42j. Verwer et al., [39] discretized in space using a basic finite difference 
approximation and employed a conditionally stable, explicit second-order Runge-Kutta-Chebyshev 
method for the time-stepping. Ang [31 formulated a dual-reciprocity boundary element method 
which is claimed to be applicable to problems in domains of arbitrary shape. More recently, in [36], 
a method involving a spatial discretization based on radial basis functions and a linearized Crank- 
Nicolson-type method for the time-stepping is formulated, while in |27l [28] differential quadrature 
is used in space, and in time a fourth-order accurate Runge-Kutta method of Pike and Roe |31] . 

Example 1. To demonstrate the optimal accuracy of the ADI scheme in various norms, we consider 
a test problem in which 

O = (0, 1) X (0, 1), Dx = D2 = \, A = \, B = 0.5, T = 1, 
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and the exact solution of ( 1.1 )-( |1.3| ) has components 

(4.1) ui{x, y, t) = cos(t) cos(27rj;) cos(7r?/), n2(x, y, t) = cos(t) cos(7rx) cos(27ry). 



Clearly g°(x,y) = [ni(x, y, 0), ^2(2^5 0)]"^, and we construct the reaction kinetic functions in (1.1) 
to be 

f(w) = f(x,y,i) + f(w), 

for any w = [wi,W2Y , where 

f(x,y,t) = ^-PAu-f(u). 



In the OSC discretization, we consider splines of degree r, where r = 3, 4, 5, and choose the 
time step r = A; = 0, 1, since it is expected that the method is 0(r^ + h'^^^~^) accurate 

in the H^{Q,) norm. For various values of N and T = 1, Tables 1 and 2 show, respectively, the 
L^(ri) and H^(Q,) norms of the errors in 1 and n/j^2; denoted by e^^i and €^^2- These norms of 
the errors are calculated using (r + 2)-point Gauss quadrature so that the quadrature error does 
not affect the accuracy of the scheme. The tables also show the experimental convergence rate of 
the scheme, denoted by Rate, calculated using the formula 



Rate 



log (error / error/^j ) 

log(/li//l2) 



where error/j 



The results in Tables 1 and 2 confirm the expected 



convergence rates of the ADI scheme in the H^{Q) norm. A; = 0, 1. 



r 


N 


l|e'i,illL2(n) 


ll<^h,2llL2(n) 


Rate 


3 


10 


0.683-04 


0.706-04 






15 


0.134-04 


0.139-04 


4.014 




20 


0.424-05 


0.438-05 


4.007 


4 


9 


0.139-04 


0.144-04 






16 


0.785-06 


0.816-06 


4.993 




25 


0.845-07 


0.878-07 


4.996 


5 


10 


0.818-06 


0.850-06 






15 


0.718-07 


0.747-07 


6.000 




20 


0.128-07 


0.133-07 


6.000 



Table 1. Example 1: L'^{n) errors at T = 1 



r 


N 




ll<^h,2llHl(n) 


Rate 


3 


9 


0.835-02 


0.864-02 






16 


0.151-02 


0.156-02 


2.975 




25 


0.399-03 


0.411-03 


2.986 


4 


10 


0.586-03 


0.609-03 






15 


0.116-03 


0.120-03 


3.999 




20 


0.366-04 


0.381-04 


3.999 


5 


9 


0.983-04 


0.102-03 






16 


0.553-05 


0.576-05 


5.000 




25 


0.594-06 


0.618-06 


5.000 



Table 2. Example 1: H^{n) errors at T = 1 
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r 


N 


ll^h,lllL°°(f2) 


\\^h,2\\ L°° (SI) 


Rate 


3 


10 


0.239-03 


0.244-03 






15 


0.472-04 


0.478-04 


4.019 




20 


0.149-04 


0.152-04 


3.985 


4 


9 


0.282-04 


0.297-04 






16 


0.162-05 


0.168-05 


4.999 




25 


0.174-06 


0.180-06 


4.997 


5 


10 


0.172-05 


0.177-05 






15 


0.152-06 


0.156-06 


5.997 




20 


0.270-07 


0.279-07 


5.983 



Table 3. Example 1: L'^{n) errors at T = 1 



With r = /i^^"*"^^/^, we also compute the L°°{Q) norm of the errors eh^i and eh,2 at T = 1 and 
the corresponding convergence rates. These L°°(Q) norms are calculated by taking the maximum 
absolute values of eh,i and eh,2 at 10 x 10 equally spaced points in each cell x [yj_i,yj], 

^ 1^ i,j 1^ In Table 3, the L°°(il) convergence rate is seen to be optimal, that is, 0(r^ + h'^'^^). 
Note that, when computing Rate in this case, we take error/^ = max{||e/i^i||Loo, ||e/i_2||L°°}- 

Invariably, OSC schemes exhibit superconvergence phenomena; specifically, one obtains an ac- 
curacy of 0(t^ + h?'^~'^) in the solution and the first spatial partial derivatives at the nodal points 
{(xj, yj)}f^=o- To examine the superconvergence of the ADI scheme, we choose r = /i*""^. From the 
results presented in Table 4, it is evident that the ADI scheme possesses the expected properties. 







Maximum nodal errors 




Maximum nodal errors 




Maximum nodal errors 




r 


N 




e/i,2 


Rate 






Rate 




deh.2/dy 


Rate 


3 


10 


0.239-03 


0.244-03 




0.951-03 


0.752-03 




0.729-03 


0.983-03 






15 


0.472-04 


0.478-04 


4.019 


0.196-03 


0.146-03 


3.897 


0.143-03 


0.203-03 


3.896 




20 


0.149-04 


0.152-04 


3.985 


0.623-04 


0.468-04 


3.982 


0.454-04 


0.644-04 


3.985 


4 


10 


0.160-05 


0.165-05 




0.969-05 


0.516-05 




0.489-05 


0.101-04 






15 


0.141-06 


0.144-06 


6.024 


0.887-06 


0.443-06 


5.895 


0.427-06 


0.928-06 


5.894 




20 


0.250-07 


0.258-07 


5.966 


0.159-06 


0.806-07 


5.979 


0.764-07 


0.166-06 


5.984 


5 


10 


0.167-07 


0.172-07 




0.986-07 


0.542-07 




0.513-07 


0.103-06 






15 


0.650-09 


0.665-09 


8.026 


0.401-08 


0.206-08 


7.896 


0.199-08 


0.417-08 


7.897 




20 


0.666-10 


0.685-10 


7.901 


0.413-09 


0.215-09 


7.906 


0.205-09 


0.428-09 


7.913 



Table 4. Example 1: Maximum nodal errors at T = 1 



Example 2. In this test problem from |38| . which has no known closed-form solution, we have 
J7 = (0, 1) X (0, 1), L>i = L»2 = 0.002, A=l, B = 2, T = 5, 

and 

g°(x,y) = [2 + 0.25y, 1 + 0.82;]'^. 

Figure 1 shows the initial (t = 0) and final (t = T = 5) graphs of u^^i and Uh^2 obtained using the 
ADI scheme with = 10, (that is, h = 0.1), r = 3 and r = h'^. The results for T = 5 (and for 
T = 10 which are similar, but not presented here) are in agreement with the theory of |38j which 
predicts that the solution u converges to the fixed point {B,A/B) when 

(4.2) 1- A + B'^ >0. 

The graphs in Figure 1 should compare with Figures 7 and 8 in ^38j and Figures 4 and 5 in 
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[36]. However, the figures in these papers are in error since they present graphs of approximations 
to ui{y,x) and U2{y,x) instead oiui{x,y) and U2{x,y), respectivelyj^j^ 




Figure 1. Example 2: The initial {t = 0) and final (T = 5) approximate solutions, Uh^i and Uh,2, 

Example 3. This test problem, presented in [T| with no numerical results, is similar to Example 
2 but with A = 3.4, B = 1, T = 1. With = 10, r = 3 and r = /i^ as in Example 2, the ADI 
approximate solutions Uh^i and u/j 2 at T = 1 are shown in Figure 2. It should be noted that they do 
not exhibit the oscillatory behavior present in |38i Figures 9, 10]. This test problem is also solved 
successfully in |36], and if Figure 8 of [36] were correctly drawn, it would be similar to Figure 2 
of the present paper. Moreover, the problem is solved in ^2] but little evidence is presented to 
demonstrate the efficacy of the method described therein. 



^E. H. TwizcU, Personal communication, July 2011. 
^ Siraj-ul-Islam, Personal communication, July 2011. 
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Figure 2. Example 3: The approximate solutions u/^^i and Uh^2 at T = 1 



With this choice of A and B, condition (4.2) is violated and, not surprisingly, convergence to 



a fixed point is not observed in numerical experiments. For example, graphs of the approximate 
solutions Uh^i and Uh,2 at T = 20 and 40 are given in Figure 3; cf., [Ml Figure 9]. To better 
understand the nature of the solutions as t increases, we plot in Figure 4 the profiles of Uh,i and 
Uh^2 at some points in O as t ranges from to 40. From this figure, it is obvious that the solutions 
are oscillatory; see also |36l Figure 9]. 







Figure 3. Example 3: The approximate solutions i and u/j 2 at T = 20, 40 
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Figure 4. Example 3: The approximate solutions Uh,i and Uh,2 against t at (0.25,0.25), 

(0.25,0.75), (0.75,0.25), (0^75,0.75)' 

Additional test problems appear in the literature. Aug [3] considered a problem similar to 
Example 2 but with 

y) = l^x'^ - -y2 - -y^f. 

However, few numerical results are presented in that paper. Verwer et al., |3,9j considered a test 
problem from |19j which has the same selection of parameter values as in Example 3 but with 

g°(x,y) = [2 + 0.25y,l + 0.8a;]^. 

They presented graphs of their approximation to U2 at various time levels. In \27\ [28], which are 
essentially identical papers, the authors claim to consider the same test problem but choose A = 1 
and B = 3.4 instead of A = 3.4 and B = 1. They also claim that graphs of their approximation 
to ui are similar to those in [39j when, in fact, [39] contains graphs of approximations to U2 only. 
Ang's test problem and Example 2 of the present paper are also considered in [271 128j . However, 
all of the graphs presented in [27\ I28| are incorrect, since, in each case, approximations to ui{y,x) 
and U2{y,x) instead of ui{x,y) and U2{x,y), respectively, are plotted. 

4.2 Gray-Scott Model 

The Gray-Scott model [HI IB] is given by ((ri])-([r3]) with 



f (u) = [F(l - ui) - ulu2, ulu2 - (F + k)u2f 
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where F and k are constants. 



Example 4. This test problem is Example 2 in [44J in which Vt 
the exact solution is chosen to be 



(-1,1) X (-1,1), r = 1, and 



ui{x, y, t) = cos(2t) cos(27ra;) cos(7r?/), U2{x, y, t) = cos(2t) cos(7rx) cos(27ry). 
Two parameter choices are considered: 

(4.3) Di = D2 = 0.001, F =1, k = 0, 

and 
(4.4) 



1, 



1, 



A; = 0. 



Zhang et al., [HI used a finite element Galerkin method for the spatial discretization with spaces 
of piecewise linear or piecewise quadratic elements defined on uniform triangulations of Q, 
and a time-stepping procedure based on a linearized second-order backward differentiation formula 
(BDF). These methods yield approximations which are second- and third-order in space, respec- 
tively, and second-order in time. At each time step, the linear algebraic system is solved using an 
iterative method based on GMRES. 



Results obtained for (4.3 ) using the ADI method with a choice of r consistent with the expected 
spatial accuracy are presented in Tables 5-7. These tables demonstrate the optimal convergence 
rates in the H^{Q),k = 0,1, and L°°(0) norms, respectively. Note that A'^, the number of space 
intervals, is chosen so that Ay, the number of time steps, is an integer. Superconvergence at the 
nodes is exhibited in Table 8 for (4.4) which is the parameter choice used in to examine rates 
of convergence. 



r 


N 


l|eh,illL2(n) 


ll<^h,2llL2(n) 


Rate 


3 


20 


0.708-04 


0.647-04 






26 


0.252-04 


0.231-04 


3.928 




32 


0.111-04 


0.102-04 


3.963 


4 


8 


0.569-03 


0.458-03 






18 


0.984-05 


0.788-05 


5.005 




32 


0.567-06 


0.459-06 


4.953 


5 


20 


0.575-06 


0.456-06 






26 


0.119-06 


0.945-07 


6.000 




32 


0.343-07 


0.272-07 


6.000 



Table 5. Example 4: L^Q) errors at T = 1 for (4.3) 



In Figure 5, we present graphs of Uh,i and u/1,2 at T = 1 computed with A = 20, (that is, 
h = 0.1), r = and r = 4. When compared with corresponding graphs (Figures 5a, 6a and 
Figures 5b, 6b) in [U], we see that the graphs are similar, taking into account their differing 
orientations. Figure 6 shows graphs of the corresponding nodal errors e^^i and at T = 1. It 
should be noted that the nodal errors in this figure are approximately 10~^ which are of the same 
order as those obtained in [HI Figures 5d, 6d]. 
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r 


N 


l|e;i,ill-ffi(n) 


l|eh,2ll/ji(n) 


Rate 


3 


8 


0.820-01 


0.722-01 






18 


0.756-02 


0.676-02 


2.932 




32 


0.137-02 


0.123-02 


2.968 


4 


20 


0.437-03 


0.358-03 






26 


0.153-03 


0.126-03 


3.989 




32 


0.669-04 


0.551-04 


3.984 


5 


8 


0.424-02 


0.346-02 






18 


0.719-04 


0.583-04 


5.029 




32 


0.404-05 


0.327-05 


5.005 



Table 6. Example 4: H^{n) errors at T = 1 for (4.3) 



r 


N 


l|eh,illL°o(n) 


l|eh,2llL°°(n) 


Rate 


3 


20 


0.175-03 


0.142-03 






26 


0.621-04 


0.506-04 


3.951 




32 


0.272-04 


0.222-04 


3.973 


4 


8 


0.106-02 


0.751-03 






18 


0.176-04 


0.120-04 


5.055 




32 


0.963-06 


0.659-06 


5.048 


5 


20 


0.102-05 


0.698-06 






26 


0.214-06 


0.145-06 


5.977 




32 


0.616-07 


0.420-07 


5.984 



Table 7. Example 4: L°°(17) errors at T = 1 for (4.3) 




Figure 5. Example 4: The approximate solutions u^^i and Uh^2 at T = 1 
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Maximum nodal errors 




Maximum nodal errors 




Maximum nodal errors 






N 












i v,aiit; 


d£h,i/dy 


deh,2/dy 




Q 
O 


90 


n OQ9 m 


n 9Q/1 OQ 
U.Z04— Uo 




n 9nQ 09 
U.zUo— Uz 


U.oy4— Uo 




U.oyD- Uo 


O 90/1 09 
U.ZU'-t— UZ 






24 


0.136-03 


0.137-03 




0.103-02 


0.432-03 


O. i zu 


0.432-03 


0.104-02 


797 
o. / z / 




28 


0.734-04 


0.739-04 


3.999 


0.543-03 


0.233-03 


4.160 


0.233-03 


0.546-03 


4.159 




32 


0.430-04 


0.433-04 


3.999 


0.327-03 


0.137-03 


3.812 


0.137-03 


0.328-03 


3.816 




36 


0.269-04 


0.270-04 


3.999 


0.201-03 


0.853-04 


4.126 


0.855-04 


0.202-03 


4.125 




90 
zu 


U.O40— uo 


U.04D— Uo 




O 90C; 0/1 
U.ZUO— U4 


U. lUo— U4 




n 1 no n/i 
U.iuy— U4 


O 90fi O/l 
U.ZUD— U4 






24 


0.115-05 


0.116-05 


fi ono 


0.721-05 


0.363-05 


^ 79S 
O. 1 zo 


0.364-05 


0.724-05 


^ 79Q 
o. / zy 




28 


0.458-06 


0.460-06 


6.000 


0.279-05 


0.144-05 


6.162 


0.144-05 


0.280-05 


6.162 




32 


0.205-06 


0.206-06 


6.000 


0.128-05 


0.647-06 


5.813 


0.647-06 


0.129-05 


5.813 




36 


0.101-06 


0.102-06 


6.000 


0.624-06 


0.319-06 


6.128 


0.319-06 


0.626-06 


6.127 


5 


20 


0.340-07 


0.342-07 




0.202-06 


0.107-06 




0.107-06 


0.203-06 






24 


0.792-08 


0.795-08 


8.000 


0.495-07 


0.248-07 


7.728 


0.249-07 


0.497-07 


7.729 




28 


0.231-08 


0.232-08 


8.000 


0.141-07 


0.724-08 


8.162 


0.724-08 


0.141-07 


8.161 




32 


0.792-09 


0.796-09 


8.000 


0.495-08 


0.247-08 


7.813 


0.249-08 


0.497-08 


7.813 




36 


0.313-09 


0.315-09 


7.873 


0.193-08 


0.984-09 


8.000 


0.985-09 


0.194-08 


8.000 



Table 8. Example 4: Maximum nodal errors at T = 1 for (4.4) 





-1 -1 



Figure 6. Example 4: The errors Ch^i and eh,2 at T = 1 



4.3 Gierer-Meinhardt Model 



The Gierer-Meinhardt model comprises (1.1)-(1.3) with 

(4.5) f(u) = [ul/u2 - ui,u\/{en) - U2/fif. 

First we examine the performance of the ADI method on test problems considered in [26] in which 
O = (—1, 1) X (—1, 1), > 0, and Di = e^, D2 = k/ fj,, where e = 0.04, fi = 0.1, and k is varied. 
Also, u(x,y,0) = [g^{x,y),g^{x,y)]'^, where 



(4.6) 5? 



1 



20 



1 + 0.001^ cos 



k=l 



kiry 



sech 



2e 



cosh ( 1 — \/x^~-\-y^] 
92 = o TTT^ 



3 cosh(l) 



In [26], this problem, which has no known closed-form solution, is solved using a Chebyshev spectral 
collocation method for the spatial discretization and a linearized backward Euler method for the 
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time-stepping. At each time step, the collocation equations are solved using a preconditioned 
GMRES method at a cost of 0(AA^/^) operations, where J\f is the number of unknowns. As in [26], 
we restrict our attention to the dynamics of ui, the initial profile of which is shown in Figure 7. 



□ 6 




Figure 7. Initial profile of 



Example 5. In this example, k = 0.0128. For r = 3, = 20 (that is, h = 0.1), r = /i^, we graph Uh,i 
at the values of t selected in [26]. Figure 8 provides graphs corresponding to Figures 7-10 in [26] . It 
is clear from our graphs that the dynamics of Uh,i follow the same general pattern displayed in [26]. 
The center of the spike at the origin sinks to the floor producing a ring formation that moves outward 
to the boundary (t = 320). Then it collapses into smaller spikes at the boundary and generates 
four new spikes in the interior near the center of each of the four boundaries {t = 420). Thereafter 
it appears that the spikes evolve to a symmetric pattern (t = 900). This pattern is more obvious in 
Figure 9 which shows aerial views of the graphs in Figure 8. There are, however, some inexplicable 
differences between our graphs and those in |'2'B] . For instance, if one compares corresponding graphs 
at t = 900, our graph shows an arrangement with four spikes along each boundary and four in the 
middle whereas in [26] there appear to be five spikes along the y-boundaries, four along the x- 
boundaries and four in the middle. Also, at t = 320, 340, the spikes near the boundary appear to 
be more prominent than in our corresponding graphs. 

Example 6. In this example, k = 0.0152 and we compare results obtained using the ADI method 
with the graphs in Figures 11-14 of |26j. In our computations, we use the same parameters as in 
Example 5, with h = 0.1, t = and r = 3. Figure 10 shows graphs of Uh^i for various values of 
t. Once again, in our graphs, we observe the general pattern obtained for Uh^i in Figures 11-14 of 
|26j . The spike splits into two spikes spreading in the x direction and becomes symmetric {t = 140). 
Then, each of the spikes splits, spreading in the y direction, and maintains symmetry [t = 290). 
Next, each of the four spikes splits into two along the x direction, and the eight spikes arrange 
themselves symmetrically about the center [t = 620). Finally, each of the four outermost spikes 
split into two and the 12 spikes arrange themselves symmetrically about the center ending with 
four spikes at the corners and eight spikes in a circular pattern around the center {t = 990). This 
is seen more clearly in Figure 11 which gives aerial views of the graphs of Figure 10. Comparing 
our pattern with that of [26j, we see a similarity up to t = 520. Kt t = 570, it appears that in 
|26j the inner spikes split as opposed to the outer ones in our case. Thus, while at t = 990 there 
are 12 spikes in our graph and in [26], the symmetric pattern that we obtain is not seen in the 
corresponding graph of [26j . 
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15 




Figure 9. Example 5: Aerial views of Uh^i at various values of t 



16 



t=30 



t=80 



t=go 




Figure 10. Example 6: Graphs of Uh^i at various values of t 
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Figure 11. Example 6: Aerial views of Uh i at various values of t 

On comparing Examples 5 and 6, we see that a modest change in the value of k leads to a significant 
change in the behavior of the solution. 

Examples 5 and 6, taken from [26J, are based on examples in |33| Section 4.2]. While the solution 
in [26] is a simple scaling of the solution in |33| (4.2)-(4.6)], the initial value prescribed is not the 
scaled initial value used in ; each of the components in the initial condition should be divided 
by e. Moreover, the graphs obtained in [26] are not given at the same t values as in [33]. In view of 
these discrepancies, the comparison given in [26j between the results presented therein and those 
of [33] is questionable. Consequently, we next consider the examples presented in [331 Section 4.2], 
where the technique employed is based on a moving mesh finite element Galerkin method with 
piecewise linear functions on triangulations of 0, with a second-order Runge-Kutta scheme for the 



time-stepping. The test problem in this case is (1.1) with f(u) given by (4.5) in which e is replaced 



18 



by e^. The initial value (which is incorrectly stated in [33^ (4.1)]; should be e) is given by (4.6). 
All other parameter choices are the same as in Examples 5 and 6. 





Figure 13. Example 7: Aerial views of u/j^i at various values of t 
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Example 7. In this example, k = 0.0128. For r = 3, N = 20 (that is, h = 0.1), r = /i^, we 
graph M/j 1 at the values of t selected in j33j. Figures 12 and 13 show the graphs and aerial views, 
respectively, corresponding to Figures 9-10 in [33J . There is fairly good agreement between the two 
sets of graphs. 

Example 8. In this example, k = 0.0152 and we compare results obtained using the ADI method 
with the graphs in Figures 11-12 of [33j. With h = 0.1, t = h? and r = 3, Figure 14 shows graphs 
of Uh^i for various values of t. From the aerial view of the solution shown in Figure 15, it appears 
that we have good agreement up to approximately t = 100. Thereafter, there is a difference in the 
arrangement of the spikes. Our figures indicate that a steady-state pattern is reached by t = 500 
which is not the case in |33j . 

Again, we see that a small change in n leads to a significant change in the behavior of the 
solution. 



t=o 



t=50 



t=100 



0.8 



V -1 -1 



X 10 



V -1 -1 



X 10 




y -1 -1 



t=500 



t=800 



t=2000 



X 10 




X 10 




X 10 




y -1 -1 



y -1 -1 



y -1 -1 



Figure 14. Example 8: Graphs of u^^i at various values of t 
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1 0.5 





Figure 15. Example 8: Aerial views of u^^i at various values of t 



4.4 Schnakenberg Model 



The Schnakenberg model j35j is given by (1.1)-(1.3) with 

f (u) = [7(a -ui+ ulu2),'y{b - ulu2)f, 

where 7,0 and b are constants. In [25], this model is studied on fixed and growing domains using 
methods based on a piecewise linear finite element Galerkin discretization in space coupled with a 
linearized backward Euler method or a linearized second-order BDF for the time discretization. The 
algebraic systems arising at each time step are solved iteratively using a preconditioned GMRES 
method. 

In the following test problems, we consider the fixed domain ^1 = (0,1) x (0,1), and take 
Di = 1, D2 = 10. 



r 


Ni 


lkh,lllL2(n) 


lk/i,2llL2(n) 


Rate 


3 


10 
15 
20 


0.683-04 
0.134-04 
0.424-05 


0.805-03 
0.159-03 
0.504-04 


3.998 
4.000 


4 


9 
16 
25 


0.132-04 
0.746-06 
0.802-07 


0.141-03 
0.792-05 
0.851-06 


4.999 
4.999 


5 


10 
15 
20 


0.775-06 
0.680-07 
0.121-07 


0.831-05 
0.730-06 
0.130-06 


6.000 
6.000 



Table 10. Example 9: L'^^Q) errors at T = 1 



0.1, b = 0.9, and construct the reaction kinetic 



Example 9. In this example, we set 7 = 10, 
functions and initial functions so that the exact solution is given by (4.1). Choosing the time step 

we present in Tables 10-13 the H^{Q),k = 0,1, and 



for various values of r as in subsection 



4.1 
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r 


Ni 


I|6h,illfri(r2) 


l|eh,2llHl(r2) 


Rate 


3 


9 


0.764-02 


0.792-01 






16 


0.138-02 


0.143-01 


2.971 




25 


0.363-03 


0.377-02 


2.993 


4 


10 


0.527-03 


0.589-02 






15 


0.104-03 


0.116-02 


3.997 




20 


0.330-04 


0.369-03 


3.999 


5 


9 


0.882-04 


0.999-03 






16 


0.497-05 


0.563-04 


5.000 




25 


0.533-06 


0.604-05 


5.000 



Table 11. Example 9: iJ^(il) errors at T = 1 



r 


Ni 


l|eh,lll-L°°(n) 


l|e/i,2llL°°{n) 


Rate 


3 


10 


0.306-03 


0.174-02 






15 


0.604-04 


0.341-03 


4.021 




20 


0.191-04 


0.109-03 


3.968 


4 


9 


0.393-04 


0.283-04 






16 


0.223-05 


0.159-04 


5.003 




25 


0.241-06 


0.171-05 


4.998 


5 


10 


0.238-05 


0.167-04 






15 


0.209-06 


0.146-05 


6.009 




20 


0.373-07 


0.262-06 


5.982 



Table 12. Example 9: L°°{n) errors at T = 1 



L°°(0) norms of the errors e^^i and eh^2 and maximum nodal errors at T = 1, respectively. These 
numerical results confirm the expected optimal rates of convergence and superconvergence at the 
nodes. 

Example 10. This test problem is Example 4 of ^25j, (cf., |34l Example 5]) in which 7 = 1000, a = 
0.126779,6 = 0.792366, and 

8 

5? = 0.919145 + 0.0016 cos(27r(x + y)) +0.01 J] cos(27rja;), 

8 

gl = 0.937903 + 0.0016 cos(27r(x + y)) + 0.01 ^ cos(27rjx). 







Maximum nodal errors 




Maximum nodal errors 




Maximum nodal errors 




r 


Ni 




Eft, 2 


Rate 




deh,2/dx 


Rate 




deh.2/dy 


Rate 


3 


10 


0.306-03 


0.174-02 




0.106-02 


0.546-02 




0.699-03 


0.989-02 






15 


0.604-04 


0.341-03 


4.021 


0.209-03 


0.106-02 


4.033 


0.137-03 


0.204-02 


3.888 




20 


0.191-04 


0.109-03 


3.968 


0.661-04 


0.341-03 


3.952 


0.438-04 


0.650-03 


3.981 


4 


10 


0.227-05 


0.166-04 




0.104-04 


0.522-04 




0.519-05 


0.992-04 






15 


0.199-06 


0.145-05 


6.024 


0.913-06 


0.452-05 


6.035 


0.458-06 


0.911-05 


5.890 




20 


0.354-07 


0.260-06 


5.965 


0.163-06 


0.816-06 


5.950 


0.823-07 


0.163-05 


5.981 


5 


10 


0.233-07 


0.167-06 




0.108-06 


0.525-06 




0.529-07 


0.994-06 






15 


0.909-09 


0.645-08 


8.025 


0.420-08 


0.202-07 


8.036 


0.209-08 


0.405-07 


7.891 




20 


0.924-10 


0.654-09 


7.958 


0.429-09 


0.205-08 


7.943 


0.213-09 


0.409-08 


7.974 



Table 13. Example 9: Maximum nodal errors at T = 1 
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In this example, we take r = 5, = 20 (that is, h = 0.1), and r = . 
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1=0.025 



1=0.125 




1=0.225 





¥ 



y 



Figure 17. Example 10: Aerial views of u/j i at various values of t 

Graphs of Uh,i at various values of t are presented in Figure 16 and corresponding aerial views 
are shown in Figure 17. On comparing the graphs of Figure 17 (rotated clockwise) with those of 
Figure 8 or 9 of [25], it is clear that there is agreement between them when t > 0.2. For t < 0.2, a 
similarity is difficult to detect as the details in the graphs of [25] are obscured because of the use 
of constant threshold shading. 



Example 11. Here we consider Example 5 of [25] where 7 = 10000, a = -0.887757,6 = 2.774242, 
and 

5? = 1.886485 + 0.001 ^°«(2vr^, gO = 0.779539 + 0.001 

i=i 3=1 

In this example, we take r = 6, = 20 (that is, h = 0.05), and r = lOO/i^ so that the number of 
time steps, Nt, for the various values of t selected in [25] is an integer. 
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Graphs of Uh^i at various values of t are presented in Figure 18 and corresponding aerial views 
are shown in Figure 19. On comparing the graphs of Figure 19 (rotated clockwise) with those of 
Figure 13 in [25], we do not observe the ripples obtained in [25j . 





□ □ 



□ □ 




□ □ 



Figure 18. Example 11: Graphs of Uh i at various values of t 





y D 



y □ □ 



Figure 19. Example 11: Aerial views of i at various values of t 
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5 Concluding Remarks 



We have presented an ADI method for solving nonhnear systems of reaction-diffusion problems of 
the form ( 1.1 ) which is second-order accurate in time and of optimal accuracy in space. The method 
enjoys several features that render it more efficient than current methods for solving such problems. 
Using test problems primarily from the literature, we have examined the accuracy of the method 
and demonstrated optimal convergence rates in standard norms for the Brusselator, Gray-Scott and 
Schnakenberg models. Also, we have compared results obtained by the ADI scheme with results 
in the literature for the Brusselator, Gierer-Meinhardt and Schnakenberg models for which exact 
solutions are not known. In the process, we have identified several errors and discrepancies in the 
literature. 

The ADI method can be generalized to solve a system of nonlinear parabolic equations with more 
than two equations. It is also possible to generalize the scheme to systems of nonlinear parabolic 
problems on rectangles where the diffusion coefficients and the reaction kinetic functions depend 
on X, y, t, u, Vu as well as when the boundary conditions are of Dirichlet, Neumann or mixed type; 
cf., [1]. In view of the high accuracy that can be achieved in space, future research will include 
the development of time-stepping methods of higher order accuracy than second; cf., 0. Also, the 
treatment of growing domains (cf., [25]) and the formulation of ADI OSC methods for problems in 
more general regions will be considered; cf., [HJE^IES]- 
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